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Abstract 

A new model of metal viscoplasticity, which takes combined isotropic, kine- 
matic, and distortional hardening into account, is presented. The basic mod- 
eling assumptions are illustrated using a new two-dimensional rheological 
analogy. This demonstrative rheological model is used as a guideline for the 
construction of constitutive relations. The nonlinear kinematic hardening is 
captured using the well-known Armstrong-Frederick approach. The distor- 
tion of the yield surface is described with the help of a so-called distortional 
backstress. A distinctive feature of the model is that any smooth convex 
saturated form of the yield surface which is symmetric with respect to the 
loading direction can be captured. In particular, an arbitrary sharpening of 
the saturated yield locus in the loading direction combined with a flatten- 
ing on the opposite side can be covered. Moreover, the yield locus evolves 
smoothly and its convexity is guaranteed at each hardening stage. A strict 
proof of the thermodynamic consistency is provided. Finally, the predictive 
capabilities of the material model are verified using the experimental data 
for a very high work hardening annealed aluminum alloy 1100 Al. 

Keywords: viscoplasticity, yield function, kinematic hardening, 
distortional hardening, rheology 
2000 MSC: 74C10, 74S05 



Nomenclature 

a distortion parameter, cf. ( fT3l) . ( 1261) 

Xk, Xd backstress and distortional backstress, respectively 

R isotropic hardening, cf. f lT8|) i 

cr cS effective stress, cf. (]20|) x 
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-RcfT 


effective radial direction, cf. (152"|L 






inelastic and elastic strains, respectively, cf. (1T1]) 1 






dissipative parts of e- u cf. ([ML, (TT4|)o 




^•ke; ^-de 


conservative parts of cf. ffT^|) 2 . ffT4jU 




A; 


inelastic multiplier, cf. fjiH|L 




V 


accumulated inelastic arc-length (Odqvist parameter), cf. (Tj 


H) 3 


s, s d 


internal variables of isotropic hardening, cf. dH]) 4 , f[2"g|) 




e 


angle between cr^ and Xd, cf. (120|) 2 




f 


overstress, cf. (TIU]). pi) 




f 


non-dimensional overstress, cf. (J6]) 




E\(K(.,a)) 


non-dimensional elastic domain, cf. ([2]) and Fig. [3b 




K(9,a) 


non-dimensional yield stress function, cf. §5§ 




R sa -\e) 


non-dimensional saturated yield stress function 




V{y,A) 


distance between y e R 2 and i C I 2 , cf. fl3]) 





1. Introduction 



Numerical simulation of complex metal forming operations is a pow- 
erful tool to reduce the development costs and to optimize the mechani- 
cal properties of a workpiece. For many polycrystalline metals, the initial 
yield surface can be approximated with sufficient accuracy by the conven- 
tional Huber-Mises yield condition which implies initial plastic isotropy. On 
the other hand, already very small plastic deformations may lead to a sig- 



nificant change of the yield surfac e compared to the initial state (lAnnin 



197 81: IWegener and Schleeel I . Il996l : iDannemever I . Il999l ; ISteck et al. 1 . 12001 



Khan et al 1 . 120 IQaJ Jbh . It is well known that the residual stresses, spring- 



back, damage evolution, and failure are highly dependent on the accumulated 
plastic anisotropy of the material. In this work we concentrate on the phe- 
nomenological modeling of the plastic anisotropy with especial emphasis on 
the distortional hardening. A conventional approach to metal plasticity is 
used in the current study: we suppose that a unique yield surface exists 
and that the material behavior is purely elastic for stresses within the yield 
surface Q 



1 Alternatively, different unconventional concepts with nu merous types of "y ield" sur 



faces exist. The so-called "subloading surface models" (see iHashiguchi I (|l989h and ref- 
erences therein) allow to capture the plastic flow for stresses within the "yield" surface. 
Such approach allows to obtain a smooth transition from the elastic into the elasto-plastic 



2 



The state of the art phenomenological plasticity is a result of accumulated 
efforts made by generations of researcher s. Unfo r tunat ely, little academic 
credit was given to the paper written by IPrager I (119351 ). Already in 1935, 
Prager combined the isotropic hardening of Odqvist type, the distortional 
hardening for the prediction of cross hardening effect, and the kinematic 
hardening for the Bauschinger effect. Interestingly, the idea of modeling 
the Bauschinger effect by the kinematic translation of the Huber-Mises yield 
surface in the stress space was taken by Prager from a conference talk given 
by A. ReuB a year before - in 1934! 



Wit hin the classical phenomenological model of lChaboche and Rousselier 
(1l983al lbh. the isotropic expansion and kinematic translation of the yield sur- 
face are considered, such that the yield surface is represented by a hyper- 
sphere in the deviatoric stress spaced Thus, the change of the form of the 
yield surface is neglected. Such models can be used to simulate the stress re- 
sponse under proportional loading. However, in general, the distortion of the 
yield surface has to be considered under nonproportional loading with abrupt 
change of the loading path. Such loading conditions are typical not only for 
multi-stage forming processes, but even for some single-stamping forming op- 
erations. In order to control the rotation of a h yperellipsoid which rep r esents 
the yield surface withing a Hill- type theory^], Baltov and Sawczuk ( 1964 ) 
introduced a polynomial representation of the corresponding 4th r ank Hill 



type anisotropy tensor in terms of the strain tensor. According to iBetten 
(119761 ). 4th and 6th rank hardening tensors are postulated as functions 



of the plastic strain. iDafalias I ( 119791 ) considered a general representation 
of the 4th ra nk tensor as a polynomia l function of the plastic strain. In 



the paper by iHelling and Miller I ( 119871 ). the 4th rank Hill- type anisotropy 
tensor is assumed to be a function of two backstress-like tensors. In con- 
trast to the above mentioned approaches, the approach of Helling allows to 



pose, 


Rccs 


(198- 


); 


Streilein ( 


1997): 


Kowalskv ( 


1999): 


Steck et al. ( 


2001 


); 


Dafalias et. al. 


2002); 


Feieenbaum and Dafalias (2007 


); 


Noman et. al. ( 


2010) 



range. Moreover, some models with a smooth stres s response can be co n struc ted using 
the concept of "bounding surface" (see, for instance, IDafalias and Popov 

2 Such yield surfa ce can be represented by a hypersphere in Ilyushin's space, as well 
(see lllvushinl (|l954 V 

3 The original approach of iHill I (|1948l ) can be used to describe a certain initial plastic 
anisotropy, but not its evolution. 
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Dafalias and Feigenbaum I (120111 ); iPietryga et. al. I (120121 ) and others mod- 



ified the Chaboche-Rousselier model introducing ordinary differential equa- 
tions which describe the evolution of tensor- valued internal variables of higher 
order (typically 4t h and 6th rank tensors) . An altern ative integral approach 



was p resented by iDanilov I (Il97ll ). In the paper of iGrewolls and KreiBig 
(120011 ). evolution equations for higher order tensors were formulated in in- 
tegral form using the Danilov's approach. Both differential and integral ap- 
proaches menti oned above allow to take the dependen ce of the hardening on 
the strain path. iKurtyka and Zyczkowski I ( 119851 . Il996j ) proposed a geometric 
approach in order to simulate a complex distortion of the yield surface. 

The rigorous proof of convexity of t he yield surface may become rather 
difficult, if the 4th rank tensors are used ( Plesek et. al. . 201 ol ). For instance , 



due to the complexity of the model presented by IPietryga et. al. I ( 120121 ). 



the convexity of the yield surface was tested numerically by plotting its two- 
dimensional projection at different loading stages. 

Probably, the most simple generalizations of the Chaboche-Rousselier 
model are based on the use of second-rank backstress-like tensors. Within 
this approach, the orientation of the yield surface follows the loading path 
such that the change of the loading direction leads to a reorientation of 
the yield surface with a certain time lag. A short overview concerning dif- 
ferent approaches based on t he use of backstress-like tensors is presented by 



pr 

Wegener and Schlegel I ( 1l996l ). In particular, within the model of lOrtiz and Popov 



( 119831 ). the size of the elastic domain along a radial line emanating from 
the origin of the yield surface depends on the angle 9 between the effective 
deviatoric stress and the backstress direction. More precisely, the critical 
norm of the effective deviatoric stress is given by a Fourier cosine series of 
9. Thus, an arbitrary yield surface which is symmetric with respect to the 
backstress direction can be approximated. On the other hand, the convexity 
of the yield surface imposes constraints on the Fourier coefficients. These 
constraints complicate the construction of practical material models, espe- 
cially if the smooth evol ution of the yield surface is intended. Another special 
case was considered by iFrancois I (120011 ). Within this approach, certain egg- 
shaped yield surfaces can be modeled, such that the egg-axis is oriented 
along a backstress-like tensor Xd and the degree of distortion is proportional 
to ||Xd||. In particular, if Armstrong- Frederick type of hardening is used to 
describe the evolution of Xd, the distortion evolves in time smoothly. The 
thermodynamic consistenc y was num e rically tested by Francois. Next, within 



the approach presented by iPanhans I (120061 ) as well as iPanhans and KreiBig 
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(120061 ). the distortional hardening was captured with the help of a tensor- 
valued internal variable of the 2nd rank. The form of the yield surface is 
given by the so-called limagon of Pascal. Within the approach of Panhans 
it can be easily guaranteed th at the elastic domain is simply connected and 
convex. Later, in the paper of IShutov et al. I ( 1201 ll ). a two-dimensional rheo- 



logical model of distortional hardening was suggested, which implies the yield 
surface to be the the limagon of Pascal. This rheological model was used to 
construct thermodynamically consistent c onstitutive equations of finit e strai n 
plasticity/viscoplasticity. In the paper by IFeigenbaum and Dafalias I ( 120081 ) . 
an existing material model was simplified in a thermodynamically consistent 
manner such that second-rank backstress-like tensors are used only. 

A relatively new concept of representative directions (see, for example, 
Freund et al. ( 201ll )) allows to generalize a uniaxial material model to cover 
an arbitrary triaxial loading. In order to compute the stress response, a 
numerical integration on the sphere S 2 is required. This concept, if com- 
bined with a uniaxial phenomenological model of plasticity/viscoplasticity, 
can produce a new phenomenological model with some realistic distortional 
hardening effects. An interesting simplified approach to the description of 



plastic anisotropy was proposed by lBarlat et. al. I (120111 ). Interestingly, this 



approach does not include the concept of kinematic hardening explicitly, but 
some distortional effects can be captured. The simplified approach to dis- 



tortional hardening, which was developed by lAretz I ( 120081 ). does not include 
kinematic hardening as well. Further, we note that some models of crys- 
tal/polycrystal plasti city allow the descr i ption of the yield surface distortion 



in a natural wa y jcf. iRousselier et. al. I tool )). For instance, in the paper 



by iFang et al. I ( 1201 ll ). the impact of microstructural hardening parameters 
on the form of the yield locus was analyzed in the finite strain context. It was 
shown that for reduced latent hardening the yield surface exhibits a larger 
curvature in the loading direction. 

A new phenomenological model of metal plasticity is proposed in the 
current study. The main features of the current model are as follows: 

(i) a two-dimensional rheological motivation of constitutive equations, which 
provides insight into main modeling assumptions; 

(ii) nonlinear isotropic hardening of Voce type and nonlinear kinematic 
hardening of Armstrong-Frederick type; 



(iii) arbitrary smooth convex yield surface for saturated distortional hard- 
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ening, which is symmetric with respect to a backstress-like tensor X^; 



(iv) degree of yield surface distortion is proportional to ||X<j||; the convexity 
of the yield surface is guaranteed at each hardening stage; 

(v) normality flow rule; pressure-insensitive plasticity; 

(vi) explicit formulation of the free energy density and thermodynamic con- 
sistency; 

(vii) overstress type of viscoplasticity according to Perzyna rule. 

In this paper, the temperature field is assumed to be constant in time 
and spacejj The model is formulated for infinitesimal strains such that 
the extreme simplicity of the current approach is not obscured by t he ge 
ometric nonlinearities. At the same time, the elegant technique of iLion 



( 120001 ). which is based on the consideration of rheological ana logies c an be 



used to generalize the constitutive equations to finite strains (IHelm 1 .12001 



Shutov and KreiBig . l2008at iHenann and Anand I . l2009t IVladimirov 1 . 120101 ) . 
As it was shown by Shutov et al. ( 201 if ), a similar technique can be imple- 
mented for two-dimensional rheological models, as well. Alternativ ely to the 



appro ach of Lion, the method of rheological models proposed by iPalmow 
( 119841 ) can be used to construct finite-strain constitutive relations. 



We conclude the introduction with a few remarks regarding notation. 
The elements of 1R 2 are denoted by x, y. The notations x ■ y := Xiyi + 
X2V2 and 



x 



Vx ■ x stand for the scalar product and the corresponding 



Itskov 


( 


2007 





Bold-faced symbols denote 



1st- and 2nd-rank tensors in Mr. Superimposed dot denotes the material 
time derivative: x = 4:X. The symbol " : " stands for the scalar product of 
two second-rank tensors 



A : B 



tr(A-B T ). 



4 The model is formulated in a thermodynamically admissible manner. Therefore, its 
generalization to thermoplasticity is straight-forward. The equation of heat conduction 
can be derived directly from the energy balance, and an additional type of free energy 
(so-called "de tached" free energy) can be i ntroduced for better prediction of temperature 
evolution, cf. IShutov and Thiemann ( 2011 ). 
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Figure 1: a) Two-dimensional rheological model is built up of a modified St.-Venant 
element (m.StV), Hooke-bodies H cxtl H^ lni i/dis, and modified Newton elements m.iVkin 
and m.Ndi s ; b) Rheological model seen from above. The angle between the (m.StV)-a,xis 
and i?dis is denoted by 9. 



This scalar product gives rise to the Frobenius norm as follows 

||A|| := VA : A. 

The identity tensor is denoted by 1. The notation A D stands for a deviatoric 
part of a tensor A D := A — |tr(A)l. 



2. Rheological analogy 

2.1. Two-dimensional rheological model 

Rheological models are useful for insight into the aspects of material 
modeling. Especially large body of information is provided b y rheological 



models if they are filled with a physical content ( jPetrov I . Il998l ). Obviously, 



the conventional 1-dimensional rheological models are not suitable for the 
description of the yield surface distortion. Therefore, all considerations of 
this section ar e carr ied out in two-dimensional space M 2 . In the paper by 



Shutov et al. I (120111 ). a two-dimensional rheological model was suggested, 
which implies that the yield surface is given by the limagon of Pascal. A new 
extended rheological model of distortional hardening will be presented in this 
section. 



In analogy to IShutov et al 



(1201 ll ). we consider a mechanical system 
which consists of a tank filled with a viscous fluid, a heavy solid which rests on 
the flat bottom (modified St.-Venant element m.StV), three elastic springs 
(Hooke-bodies H ext , H kin , and i^dis) connected to the modified St.-Venant 
element, and two spheres (modified Newton elements m.A^ in and m.Ndi S ) 
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Figure 2: Behavior of idealized two-dimensional bodies: a) Hooke-body H. Its elongation 
is given by Eh = R 2 ; b) Modified Newton-body m.N. Inelastic arc-length p is 

used instead of the physical time t to formulate the constitutive equations; c) Modified 
St.-Venant element m.StV. The friction depends on the angle 6. 



floating on the surface of the fluid (Fig. [T])J3 The mechanical properties of 
these idealized bodies are postulated as follows: 

• (H): For the Hooke-bodies (see Fig. [2k), the spring force au is pro- 
portional to the length of the body, and the force is oriented along the 
spring axis: cx# = c Eh- Here, e# = M , and c > is the stiffness 
of the spring. In particular, the Hooke-body possesses a zero length in 
unloaded state. 

• (m.N): The two-dimensional Newton element is represented by a sphere 
which is floating on the surface. Following the Newton's law of viscous 
flow, we assume that the fluid resistance to the motion of the sphere 
is proportional to its velocity -^En- Thus, 4eV — K &n, where x > 
is a viscosity parameter. Next, in order to obtain the constitutive re- 
lations of the modified Newton-element, the physical time t is formally 
replaced by the accumulated inelastic arc-length (Odqvist parameter) 
p. Thus, we postulate for (m.iV)-element (see Fig. [2b): 

^-e n = x a N . (1) 
dp 

Such modification is possible whenever the inelastic arc-length p is 
available. The arc-length p will be introduced formally in the following. 
The use of this parameter instead of the time t allows to construct r ate- 



independent constitutive relations (see, for example. lHaupt I ( 120021 )). 



5 An animated version of the rheological model with only one modified Newton element 
is available at |http://www.youtube.com/watch?v=QEPc3pixbC0 
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• im.StV): The heavy solid rests upon the bottom of the tank and there 
is a friction between them. By a, —x^, and — Xd denote now the forces 
acting on the (m.StV)-element due to the elongation of the Hooke- 
bodies H e ^ t , ifkin, and Hdi S , respectively (see Fig. [2b). The force a will 
be understood as an external load; Xk and Xd will be responsible for 
the effects similar to kinematic and distortional hardening, respectively. 
The resulting (effective) force is thus given by <r e ff = a — x k — x&. Let 
the axis of the (m.SYV^-element be always oriented along the resulting 
(effective) force a e g. We suppose that the fluid resistance opposed the 
rotation of the solid is negligible. The (m.StV)-element remains at 
rest as long as ||oW|| < y2/3K, where ^/2/3K > is a nonconstant 
friction. The function K is computed as follows. For Xd = 0, we put 
K = K , where K is a given basic friction. Further, suppose that Xd ^ 
0. Let 9 be the angle between the axis and Xd- 9 = arccosf ,,Z ca n n I ■ 
Moreover, let a = ||xd||/a;™ ax be a distortion parameter, which is a 
unique function of \\xd\\- Here, x™ ax > is the upper bound for ||xd||, 
therefore we get a G [0,1]. Finally, we consider the friction to be 
a function of 9 and a: K = K(9, a) Kq. In particular, for a fixed 
Xd 7^ 0, the friction K depends solely on the angle 9. A simple ansatz 
for K(9, a) will be presented in the next subsection. 

Remark 1. The choice of notations in this section is dictated by the need 
to keep the structure of the rheological model similar to the structure of small 
strain plasticity. For that reason, the forces imposed on the (m.StV)-element 
by the Hooke-bodies are denoted by — x^, —Xd rather than x k , x d - 

2.2. Direction- dependent friction and definition of overstress 

Let e\ = (1, 0) G M 2 . In this subsection we construct the non-dimensional 
function K(9, a) which plays a central role in the current study. In terms 
of the rheological model introduced above, this function is understood as a 
friction coefficient, but in the following sections it will be treated as a non- 
dimensional yield stress. It is useful to interpret such functions geometrically 
in terms of a parametric family of closed subsets in R 2 : For each a G [0, 1] the 
corresponding subset El(^(-, a)) consists of y G M 2 such that \\y\\ < K(9, a), 
where 9 G [0, 7r] is the angle between y and e*i (cf. Fig. [3b). More precisely, 
we put 

E\(K(-,a)):={yeR 2 /{0} : \\y\\ < K(9, a), where 9= U{0}, (2) 
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Figure 3: a) A closed unit disc which is used to represent the undistorted yield surface; b) 
Geometric interpretation of K(9, a) for a fixed a. For each a a closed set El(K(-,a)) C K 2 
is defined; c) A convex set El sat corresponds to the given function K sat (9). The boundary 
of El sat can be associated with a saturated form of the yield surface. 



(y, ei) := arccos . 

MMI ll e ill y 

First, for a = 0, we postulate that 0) = 1, which implies that El(K(-, 0)) 
is a closed unit disc in M 2 (see Fig. [3b,). Further, suppose that a smooth func- 
tion K s *\6) is given such that El sat := El(i? sat (-)) is convex and i? sat (0) = 
1 (see Fig. [3b). We need to find a smooth function K(9,a) such that 
E1(X(-, a)) will be convex for all a and 

tf(0, a) = 1 for a G [0, 1], K(6, 1) = K sat (6) for G [0, tt]. 

Remark 2. As it will be shown in the following, the given function 
K sat (6) corresponds to the form of a saturated distortional hardening with 
the maximum distortion (Fig. [3b). 

Remark 3. The parameter a should be understood as a distortion pa- 
rameter, such that a = and a = 1 correspond to zero and maximum 
distortion, respectively. 

In other words, an interpolation rule is needed between the intact ini- 
tial unit disc (corresponds to a = 0) and the maximum distorted set El sat 
(corresponds to a = 1). 

Remark 4. Unfortunately, the linear interpolation rule K(8, a) = (1 — 
a) + aK sat (8) is not suitable, since, in general, the convexity of El(K(-,a)) 
is violated for some a G (0, 1). 

The interpolation rule which is implemented in the current study is con- 
structed as follows. First, for any y G M 2 , A C M 2 we define the distance in 
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Figure 4: a) Definition of the distance T)(y,A) according to (|3|); b) Definition of y(9,a) 
for 9 £ [0, 7r], a € [0, 1]; c) A line of constant overstress / for a fixed a £ [0, 1]. 

a natural way (see Fig. 

V(y,A) := inf \\y-z\\. (3) 

Next, we define the product of a G IR and A C M 2 as := {ay : y e A}. 
The set El(^(-,a)) is obtained from the set aEl sat by adding additional 
points whose distance from aEl sat does not exceed 1 — a (see Fig. Hb) 

El(K(-,a)) := {y G R 2 : P(f,aEl sat ) < 1 -a}. (4) 

Since El sat is convex, so is El(^(-,a)). 

Formally, since El sat is convex, for each 9 e [0,7r], and a G [0,1] there 
exists a unique y(0,a) such that y = ||^T|| (cos(6 ) ), sin(6')) and V(y, aEl sat ) = 
1 — a (see Fig. Hb). Thus, in accordance with we put 

K(0, a ):=\\p{e, a )\\. (5) 

In what follows, each set El(^(-,a)) will be used to reflect the elastic 
region in the stress space. Since a viscoplasticity model of overstress type is 
to be constructed, a proper definition of the overstress will be needed. For 
given y G 1R 2 , a G [0, 1] we define a non-dimensional overstress as a distance 
from elastic domain: 

f(y,a) :=v(y,E\(K(;a))) = (V(y, aEl sat ) - (1 - a)), (6) 

where (x) := max(x,0). The definition is summarized in Fig. Hfc. Thus, 
for the numerical computation of the overstress it is sufficient to evaluate 
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A sat (7r) < 1 = A' sat (0) A sat (vr) = 1 = A sal (0) A' sat (7r) > 1 = A' sat (0) 



Figure 5: The material function K sat (9) is not uniquely determined by the form of the 
saturated yield surface. The function is unique if K aat (n) > is additionally specified. 

T>(y, aEl sat ). Such computation can be performed explicitly if the boundary 
of El sat is represented by a set of circular arcs. A concrete algorithm is 
presented in Appendix A. Note that such requirement is not restrictive, since 
any smooth convex curve can be approximated by circular arcs with sufficient 
accuracy. In most cases 4 or 5 arcs are sufficient for practical purposes. 

Remark 5. The definition fl6]) of the overstress / will be advantageous 
in connection with a normality flow rule. In particular, the derivative of the 
overstress with respect to y possesses a unit norm for / > 0, i.e., || 9 ^q^ || = 1. 

Moreover, the set {y : f(y) < fo} is convex for all fo > 0. 

Remark 6. Using the interpolation rule proposed above, the function 
K(6, a) is uniquely defined by the material function K sat (8). Here, the func- 
tion K sat (6) describes the saturated form of the convex symmetric yield sur- 
face. In some cases this form can be identified experimentally (cf. Remark 
9). Figure demonstrates that the function K sat (8) is not uniquely deter- 
mined even if the form of the saturated yield surface is known. This is due 
to the fact the the position of the origin relative to El sat is not unique. The 
function K sat (8) is uniquely determined by specifying K sat (7i) > 0, which 
is a material parameter "hidden" in the material function K s&t (6). This 
parameter should be chosen in such a way that a better description of the 
experimental data is achieved. 

2.3. Some constitutive equations in two-dimensions 

Consider a system of (constitutive) equations as follows. The total dis- 
placement of the point A in Fig. [I] with respect to the bottom will be denoted 
by e. This displacement is a sum of the elastic elongation of the (if ext )-body 
and the inelastic displacement of the (m.SW)-body, denoted by e e and ej, 
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respectively: 

e = e c + ej. (7) 

The displacement of the (m.StV)-body, in turn, is composed of the elastic 
(ifkin)-elongation and the inelastic (m.A^ in )-displacement. Analogous de- 
composition holds for Hdis and m.A^is, as well: 

£*i = ^kc + £*ki, £i = £*de + £*di- (8) 

The total potential energy of the system equals 

if) = lf>el(e e ) + V'kin(eke) + V'dis^dc) = A* ll^cf + y H^kcf + y 1 1 ^dc 1 1 2 • 

Here, /x, Ck, Cd > are the stiffnesses of H ext , H kin , and i^dis, respectively For 
the forces a, x k , and x& we get 

dip e \(e e ) _ o^wnCek) - d^dis^de) , n s 
de e de ke de de 

The overstress / is defined as a function of ||<? e fr|| = \\<r — x k — s d ||, 9 = 
arccos f-pr^p^ J , and the distortion parameter a = WxdW/x™* by 

f(a,x k} x d ) = f(\\a eS \\,6,a) := \J\ K Q f ( 10 ) 
" a< '"'" (cos(0),sin(0)). (11) 



Due to the fact that || ^Wf 1 ) || _ j f or positive overstress, we have 



df(a,x k ,x d ) 



d, 



a 



1, for / > 0. 



We postulate the normality flow rule (normality to the hypersurface of con- 
stant overstress ) in combination with the Perzyna-type of viscoplasticity 
flPerzvna I . fl~963h 

£i = Ai — , for / > 0, £i = for / = 0; A t = — f — / 

Oa T]\Ko 

Here, rj > and m > 1 are parameters of the Perzyna rule; fco > is used to 
get a dimensionless term in the parentheses. 
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Note that the forces x k and x d act on the modified Newton-elements 
m.N kin and m.N dis , respectively. Thus, we get in accordance with (jTJ 

U £ k ; -. d Edi 

— = X k X k , —. = X d X d - 

dp dp 

Here, x k , x d > are modified viscosity parameters describing m.N^ and 
■m.Ndis, respectively. Let the evolution of the inelastic arc-length p be given 
by p = A; = ||£i||. Thus, we obtain 

e ki = Ai x k f k , e di = Ai x d £ d . (12) 

It follows from (j^g and (fl2"]L that for proper initial conditions we have 
ll^dll < V x d- By putting x™ ax = l/x d we specify the definition of the 
distortion parameter a (cf. Section 2.1) 

a := x d ||x d ||. (13) 

Equations f lT2|) in combination with ()9]) 2 and OHJJg describe the evolution of 
the "backstresses" in the hardening/recovery format. The saturation of x d 
implies the saturation of the "distortional hardening", which takes place 
much faster than the saturation of the "kinematic hardening" . Thus, "slow" 
and "fast" saturation should be assumed for x k and x d , respectively. Note 
that the dependence on the strain path is captured in a vivid way, such that 
the system exhibits fading memory: only the most recent part of the e^-path 
influences the current "stress" state a. 



3. Material model of viscoplasticity 

3.1. Closed system of constitutive equations 

Let us formulate a system of constitutive equations of viscoplasticity. 
First, we suppose that the volumetric response is elastic. More precisely, 
the hydrostatic stress component trcr is assumed to be a linear function of 
trs. Next, suppose that the deviatoric stress component cr D depends solely 
on the history of the strain deviator £ D . In order to describe this depen- 
dence, we generalize the two-dimensional constitutive equations presented 
in the previous section to five dimensions^ During the generalization we 



6 Mathematically, cr and e are elements of a 5-dimensional vector space of trace-free 
(deviatoric) second rank symmetric tensors. 
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have to make sure that the resulting model inherits the properties of the 
two-dimensional rheological model. Toward that end, the displacements and 
forces are formally replaced by deviatoric strains and stresses, respectively; 
the scalar product in R 2 is replaced by the scalar product of two second-rank 
tensorsQ In order to take the isotropic hardening into account, the constant 
parameter K is now formally replaced by K + R, where R is a harden- 
ing variable. In order to describe the evolution of R, we introduce a scalar 
strain-like internal variable s (which is similar to the inelastic arc- length p), 
its dissipative part s d , and its conservative part s c . 

For the strain tensor e consider its inelastic part S\ and elastic part e e . 
Let £ki and £ ke be the dissipative and conservative parts of e h which are 
connected to the nonlinear kinematic hardening. Analogously, and e dc 
are parts of S\ associated to the distortional hardening. More precisely, we 
postulate 

e = e c + Si, Si = e kc + e k i, £i = £dc + £di, s = s e + s d . (14) 

Note that the first decomposition is related to (J7J). Moreover, (fl4"]) 2 and 
(JHJ3 can be motivated by (jHJ). The evolution of the state of the material 
is captured by the inelastic flow and the inelastic flow (s ki , e di , s d ) which 
takes place on the microstructural level. 

The specific free energy per unit mass is given by 

tj) = 1p e l(e e ) + lp khl (e ke ) + ^dis^de) + ^iso( S e), (15) 

pM?e) = |(tr s c f + /x ||6 c D f, p^ km (s ke ) = I ||s k D J 2 , (16) 

MiisOde) = y Ikdef. P^iso(Se) = T^) 2 - ( 17 ) 

Here, k, [A, c k , c d , 7 > are material parameters; p > stands for the mass 
density. The quantity ip c i(e e ) stands for the energy stored due to macroscopic 
elastic deformations. The remaining part ip kin + ?/> dis + ^ iso is used to cap- 
ture the energy associated with the defects of the crystal structured Next, 



7 Using a similar approach, a two- dimensional rheologic al model was already generalized 
to cover finite strain vsicoplasticity (jShutov et al. ~l. l201lh . 

8 Note that i/ikm + 4>dis + V'iso does not necessarily reflect the entire "defect energy" . It 
is natural to assume that a part of the "defect energy" is not connected to any hardening 
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we postulate the following relations for stresses, backstresses, and isotropic 
hardening 

dip e i(e e ) di> kin (£ ke ) dip dis (e de ) dip iao (s e ) 

cr = P — ^ , Ak = p — , Ad = p — , K = p — . 

de e de ke de de ds e 

(18) 

Substituting (PBJ and (HTJ) into ([IS]) we get 



(7 



k tr(e e )l + 2//e°, X k = c k e^ e , A~ d = c d e% e , R = 7 s e . (19) 



On the one hand, these relations can be motivated by the rheological 
model from Section 2. On the other hand, as it will be shown in the following, 
relations (TIB"]) will be sufficient for the thermodynamic consistency of the 
material model. It follows immediately from ( !T9|) that tiX k = trX d = 0. 

Suppose that the degree of distortion a depends solely on ||A~ d ||. A 
concrete dependence will be specified in the following (cf. (j26p ). The effective 
stress tensor and the angle 9 are defined now through 

<r eff := cr - X k - X d , 9 := arccos ( °" C p - ^ d ) . (20) 

V ll Cr effll ll A d|K 

Note that for \\X d \\ = the angle 9 is arbitrary. To be definite, we put 6 = 
in that case. Further, analogously to (TTU]) and (ITT]) , we define y e M 2 and 
the corresponding overstress / (see Fig. Ek) 

/!>, X t , X d , R) = /(||<rS[ll. ». «. : = yft^o + R) Ma), (21) 

?: =vpiU (cosW ' sinW) ' (22) 

Here, > is a fixed material parameter (initial yield stress), and the 
function f(y,a) is defined through (J6j). 



mechanism ( Shutov and Thiemann . 2011 ). Alternatively. iHenann and Anand" (2009) pre- 
fer to neglect the free energy storage 'iAiso which is associated with the isotropic hardening. 
Interestingly. iFeigenbaum and Dafalias I ( 2008 ) suggest that the defect energy is released 
while the distortion of the yield surface takes place. Thus, again, somewhat smaller energy 
storage will be predicted than by assumption (j 15[) . The choice between many alternatives 
is important for the prediction of the inelastic dissipation and should be based on relevant 
experimental observations. 
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Deviatoric stress space K 2 Deviatoric stress space 

Figure 6: a) Sketch of the yield surface in the deviatoric stress space and definition of the 
overstress /(<r) (cf. ([20]) , pTjl ). The elastic domain in the stress space is associated with 
El(K(-, a)) Cl 2 ; b) Sketch of the proof of the inequality (j27]) . 



The elastic domain corresponds to stress states with zero overstress /. For 
a given stress tensor cr, a non-dimensional vector y e M 2 must be evaluated 
according to (1221) . Observe that the angle between y and e*i coincides with 
the angle between cr® s and Xd (see Fig. Ek). According to (ED), the stress 
state er lies within the elastic domain if and only if y G El(K(-,a)). The 
origin of the elastic domain corresponds to {y = 0} or, equivalently, {cr® s = 
0} = {cr D = Xk + Xj}. Next, observe that the direction of the elastic 
domain coincides with the direction of X^, and the size of the elastic domain 

in that direction equals tJ~^(K + R) (see Fig. [6^)- 

Remark 7. The relations (l2Tj) and fl22l) imply that the resistance to 
plastic deformation depends on the actual direction of the loading relative 
to the recent loading path. This dependence on the loading direction can 
be associated with the activation/deactivation of crystallographic slip planes 
as well as mobilization/demobilization of oriented dislocation structures. It 
is well known that the gliding of dislocations is obstructed by the cell walls 
under monotonic loading (6 = 0), but after a strain path ch ange (0^0), the 



loadi ng may drive the disloc ations toward th e cell interior (iViatkina et al. 



20071 ). In the monograph by I Viatkina I (120051 ). the mechanism of "directional 



remobilisation" under the strain path change is explained as a remobilisa- 
tion of dislocation locks and dipoles which were formed during the previous 
loading, in contrast to the statistical remobilisation which is independent of 
the strain path change. 
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The normality flow rule in combination with the Perzyna rule is used 



. df(tr,X*,X d ,R) 1 f 1 A m 
Si = A ; for / > 0, Si = for / = 0; Ai = - -z-f . 

OtT T] V Kq J 

(23) 

Here, Ai > is an inelastic multiplier which controls the rate of the inelastic 
flow. Indeed, since || 9 \\ = 1 for / > 0, we get (cf. (H3|) ) 



9?7 

df((T,X k ,X d ,R) 



da 



lfor/>0, ||ei||=Ai, p = \i. (24) 



We emphasize that k is not a material parameter, and we put k = 1 MPa. 
A concrete algorithm for the evaluation of the derivative is presented 

in Appendix B. Note that the normality rule ( 123]) implies an incompressible 
flow: tr£i = 0. In order to take the saturation of the kinematic and dis- 
tortional hardening into account, we postulate for the inelastic flows on the 
microstructural level (cf. ( TP2"]) ) 



e k i = A; x k X k , e di = A; x d X d . (25) 

Here, x k , x d > are material parameters. Recall that tr X k = tiX d = 0. 
Thus, the inelastic flow on the microstructural level is incompressible as well: 
trSki = tr^di = 0. It can be easily shown that for U-X^ll^o — we have 
ll^dll < l/ x d- Thus, analogously to f lT3|) . we define the distortion parameter 
a E [0, 1] by 

«:=x d ||X d ||. (26) 

For a given deviatoric stress cr D consider a convex set & = {er* : trcr* = 
0, /(""*) < /( " D )}- The gradient d ^ CT ' x ^ Xd ' R ^ coincides with the unit 
outward normal n to the boundary of & at <r D . Moreover, the state cr* = 
Xk + X d lies within (3 (cf. Fig. |6b). Due to the convexity of & we have 

of ff : 6i H Ai(<r D - tr*) : n > 0. (27) 

Having this inequality in mind we formulate the evolution equations for the 
internal variable s and its dissipative part s d . 

cr eff : £; tre i= o cr° ff : £; . /3 

S = 7^ = K^H< Sd= 7 Sii • (28> 
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where f3 > is a material parameter controlling the saturation of the isotropic 
hardening. It follows from (|2"T|) that s > 0. Thus, similar to the inelastic arc- 
length, the variable s increases monotonically. Note that in the case of pro- 
portional monotonic loading we have 9 ~ 0. Thus, (t^ s : &\ ~ Hc^H ||£i|| = 
\\cr® s \\ A;. Moreover, for slow loading we get / <C (K + R). Therefore, 
llcr^H « a/2/3 (i^o + -R)- Thus, under quasistatic proportional loading, the 
parameter s evolves similar to the inelastic arc- length: s « ||ei|| = 

y/2/3\{. Under general loading conditions, although, the evolution of s de- 
pends not only on the rate of the plastic flow, but also on its direction. 

Finally, the system of constitutive equations is closed by initial conditions 
imposed on the strain-like internal variables 

_ I _ ^0 _ I _ _0 _ I _0 „| _ _0 _ I „0 

£i\t=0 — &i j EjdU=0 — £ ki; £ di|t=0 — £ dii s \t=0 — 8 , Sd |t=0 — S d . 

We suppose tref = tre^ = tre^ = 0. If the undeformed state is assumed 
to be stress free at t — 0, then e? = 0. The quantities £ ki and e di can be 
used to capture the initial plastic anisotropy of the material]^ In particular, 
the yield condition at t — does not have to coincide with the Huber-Mises 
criterium. 

3.2. Proof of thermodynamic consistency 

Let us cons i der th e Clausius-Duhem inequality in the form (see, for ex- 
ample, lHaupt I ([2002b) 

Si := -a : e - ip > 0. (29) 
P 

Taking the kinematic relations (fl4|) into account, we rewrite the stress power 
as follows 



a : e = 

a : (e e + ei) - X k : ej + X k : (e ki + e kc ) - X d : e ; + X d : (e di + e de ). (30) 
Moreover, differentiating (fT5|) . we get for the time derivative of the free energy 

• dlp e l(E e ) . 9^ kin (£ k e) . <9^ d is(ede) • , <9^iso(s e ). / 01 x 

V = — n : £ e H ?; : £ke H k : S de H ?; s e . (31) 

de c de ke de dc ds e 



This is equivalent to the introduction of initial backstresses. 
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Substituting (I3"U|) and ( T5TT) into (12T)|) and taking the potential relations ( TTS1) 
into account, we obtain the Clausius-Duhem inequality in the following form 

p 5i = (<T e $ ■ Si - R s) + X k : e ki + X d : £ di + Rs d > 0. 

It follows immediately from (1231) and (I2"gj) 2 that X k : £ ki > 0, X d : £ di > 0, 
and Rsd > 0. In order to prove the thermodynamic consistency of the 
material model it remains to show that cr c fj : e\ — R s > 0. Indeed, 

<r cff : ei - R s ^ <x c D ff : £i (1 - R/(K + i?)) = <x° : s ; (iT /(^o + R)) > 0. 

The thermodynamic consistency of the material model is thus proved. 

Remark 8. Note that the proof of the thermodynamic consistency is 
essentially based on the inequality cr® s : k\ > 0. Any flow rule which gov- 
erns £; and complies with this inequality would yield a thermodynamically 
consistent material model, as well. For instance, the radial flow rule can be 
considered as a simplified alternative to the normality rule fl23l) 1 

_D 

ti = XiR c s, R c {[ := Tj — g-77. (32) 

ll°"effll 

3.3. Identification of material parameters 

The material model contains 11 material parameters and a material func- 
tion K sat (9). Let us discuss the identification of these quantities. First, the 
elasticity parameters k and fi can be determined basing on the experimental 
data for elastic deformations. Next, the initial yield stress Kq can be cali- 
brated using the graphical method from a quasistatic uniaxial tension test. 
The viscosity parameters rj and m of the Perzyna law are typically identi- 
fied using a series of tests under monotonic loading with different loading 
rates. Further, the material function K SAt (6) is uniquely determined for the 
fixed ^ sat (7r) if the form of the saturated yield surface is known (for details 
see Remark 6). For simplicity, one may assume ^ sat (7r) = 1. In that case, 
the parameters of isotropic hardening (7 and (3) can be identified using the 
information about how the size of the elastic domain evolves under mono- 
tonic loading. Finally, it remains to identify two parameters of kinematic 
hardening (c k and x k ) and two distortional parameters (c d and x d ). This 
can be done by minimization of a least-squares functional which represents 
the discrepancy between measurements data and corresponding model pre- 
dictions. Experimental measurements related to non-proportional loading 
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are necessary in order to obtain a reliable identification procedure. Some 
regularization techniques can be used to reduce the correlation among the 
parameters and to reduce th e probability of getting trapped in local minima 



flShutov and KreiBig I . l2010f ). 

The success of the identification procedure depends on the quality of 
initial approximation chosen for the unknown parameters Ck, Cd, Xk, and 
Xd. The order of magnitude of these parameters can be estimated basing 
on the following considerations: The upper bounds for ||Xk|| and ||Xd|| are 
given by x^ 1 and x^ 1 , respectively. At the same time, the increment of 
the inelastic arc-length which corresponds to the saturation of kinematic and 
distortional hardening under proportional loading is proportional to (ckXk) -1 
and (cdX d ) _1 , respectively. 

4. Numerical computations 

In this study, for simplicity, the evolution equations (1251) . ([25]), and (1251) 
are integrated numerically using explicit time-stepping scheme. If rate- 
independent material response is to be simulated, a viscous regularization 
with fictitious small viscosity rj > can be used0 

In this section we validate the predictive capabilities of the material 



model. Toward that end, we consider experimental data of iKhan et al. 



(l2010al ) obtained for a very high work hardening aluminum alloy - annealed 
1100 Al. The yield points were identified experimentally under combined 
tension-torsion of thin-walled tubular specimens using a small proof strain. 
In order to simulate the deformation of thin-walled tubular specimen we 
compute the stress response at a single material point. Consider a Cartesian 
coordinate system such that its basis vectors ei, e2, and e3 are oriented along 
the local axial, hoop, and radial directions, respectively. The stress state can 
be idealized approximately as a special case of the plane stress: 

cr = (j n ei <g> ei + er 12 (ei <E> e 2 + e 2 <g> ei). (33) 

Here, an and ayi are associated to the axial and torsional loading, respec- 
tively. The measurement results are represented in the (an, v / 3cr 12 )-space in 
Fig. Uh- The initial yield surface can be approximated with sufficient ac- 
curacy using the conventional Huber-Mises yield condition. Therefore, an 



10 Note that such viscous regularization allows to smoothen the sharp transition between 
elastic and plastic regions. 
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initial plastic isotropy will be assumed during the material modeling. The 
yield points which were determined after 2% axial prestrain are depicted in 
Fig. [TJi, as well. The material parameters used to simulate the material 
response are summarized in Table [Tlf"1 Moreover, due to the initial isotropy 
we consider the initial conditions as follows: 

c-0 o o _ n o _ o _ n 

The smooth function K sat (9) which is needed to compute the overstress f(y) 
corresponds to the saturated form shown in Fig. [3b. Due to the axial pre- 
strain, the distortion parameter a ranges from up to 0.999993, which cor- 
responds to the (almost) saturated distortional hardening. As it is shown in 
Fig. [7Ji, the yield locus undergoes the isotropic expansion, kinematic trans- 
lation and distortion. The surfaces of constant overstress are depicted in Fig. 
[7b. In accordance with the modeling assumptions, these surfaces are slightly 
less distorted than the corresponding yield surface. 

Remark 9. Note that the form of the yield surface in the (an, Vocr^)- 
space coincides with the boundary of E\(K(-,a)). This is due to the well- 
known fact that the scalar product of two symmetric tensors cr^ and cr n 
(where crj and crjj comply with fl3"3"j) ) corresponds to the product of two 
vectors aj,an E R 2 , defined by 07 := (a[ x , y/SaL), &11 '■= (&[{■> \/3crf|). 



More precisely 

CT ? : = g^Wi + Scr^oia) = g<?/ • on- 

Therefore, in the context of ([3B1 . the angle between two deviatoric stress- 
states coincides with the angle between two corresponding vectors in the 
(an, v / 3o"i2)-space: 

11 \ f o-j ■ a n 



arccos | - — g-— - — j^ 77 ' = arccos 



o~i || 11°" //IK V ll cr /Il 11°"// 

Thus, if the form of the saturated yield surface is determined experimentally 
in the (an, v3<7i2)-space, it can be used to identify the smooth function 
K sat (9). 



n It is not the aim of the current study to identify the material parameters corresponding 
to the 1100 aluminum alloy. Instead, we validate the material model by the qualitative 
description of the real experimental data. 
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Figure 7: a) Experimental data for annealed 1100 aluminum alloy (jKhan et al. I . l2010af ) 
and corresponding simulation results, b) Lines of constant overstress / corresponding to 
plastic anisotropy induced by 2% prestrain in the axial direction. 



Table 1: Material parameters used to validate the material model: the viscosity effects are 
neglected. 
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If, additionally, the specimen can be loaded by an internal pressure, the 
hoop stress <r 22 must be considered, as well: 

cr = a n ei <g> e x + a 22 e 2 <S> e 2 + o- 12 (e! <g> e 2 + e 2 <g> ex), cr 22 > 0. 

We simulate the evolution of the yield surface in the process as follows. 
Starting from the same isotropic initial state, a 2% prestrain is prescribed 
in the hoop direction. Thus, a similar plastic anisotropy is introduced, as 
in the previous case. The form of the yield surface for a fixed axial stress 
on is then represented in the (<r 22 , v / 3o"i2)-space. As it can be seen in Fig. 
El the form and the size of the yield loci for cr n = 5 MPa, er n = 10 MPa, 
and an = 15 MPa are the same as for an = —5 MPa, an = — 10 MPa, and 
on = —15 MPa, respectively. Similar to the conventional Huber-Mises yield 
condition, the material yields at larger <r 22 stresses for positive en than for 
negative a u . 
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Figure 8: Influence of the axial stress <j\\ on the form of the yield surface after 2% prestrain 
in the hoop direction: a) Positive axial stresses; b) Negative axial stresses. 



5. Conclusion 

A new material model of metal viscoplasticity with an extremely sim- 
ple structure is presented in the current study. The main modeling as- 
sumptions are visualized with the help of a new two-dimensional rheolog- 
ical model. Only second-rank backstress-like tensors are used to capture 
the path dependent evolution of the plastic anisotropy. Thus, mathemati- 
cally, the model is not much mor e complicated than the classical model of 



Chaboche and Rousselier I fll983al Jbl). At the same time, the proposed tech- 
nique possess considerable generality and flexibility. No specific form of the 
saturated yield locus is considered, since any smooth convex yield locus can 
be captured. 

An important ingredient of the material modeling is the interpolation be- 
tween the initial intact yield surface and the saturated one. The interpolation 
rule proposed in the current study ensures the convexity of the yield surface 
at any stage of hardening. It is shown that this interpolation rule allows to 
obtain a thermodynamically consistent material model. 

The model contains 6 hardening parameters with 2 parameters per hard- 
ening type. These parameters posses a clear mechanical interpretation and 
can be identified experimentally. 
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Appendix A 

Let us discuss the numerical computation of T>(y, aEl sat ). Suppose that 
the upper half of the boundary of El sat is given by N circular arcs connecting 
y 4-1 and y 1 for i £ {1,...,N} as shown in Fig. [Hh- The outward normal 
and the tangent at y 1 will be denoted by n 1 and t % , respectively. The 
orientation of the tangent is chosen in such way that the pair {n \ t *} forms 
a right-handed corner. In particular, we have 

y° = (i,o), r° = (o,i), r^ = (o,-i). 

For each arc connecting y i_1 and y 1 consider its center y c l and its radius r\ 
In order to make sure that the boundary of El sat is smooth, we require 

->0 -» 1 -> N -» N 

y - Vc _ H o v - yJ _ n N. 

11 ' ~ - V = n i for all i £ {1, JV - 1}. 



Let a £ [0,1] and y = \\y\\ (cos(#), sin(6')) be given such that 9 £ [0, 7r]. 
The algorithm used to compute the distance T>(y, aEl sat ) is as follows. 

• Check the inclusion: If \\y\\ < aK sat (9) then £>(?/, ctEl sat ) = 0. 

• Otherwise, find the corresponding arc: Find % £ {1,...,N} such that 
(cf. Fig. Eb) 

(y-ay 1 - 1 ) ■ P' 1 > 0, (y- ay') ■ > 0, (34) 

(y-ay^-Q-iy' -y*- 1 ) > 0, Q := ei ® e 2 - e 2 ® el. (35) 

• The distance is then given by 

£>(y,aEl sat ) = \\y- ayJW - ar\ 

Moreover, the outward unit normal is given by 

H = df(y,a) = dV&aEr*) = y-ay c > 
dy dy \\y - ayj\\' 
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Figure 9: a) A smooth boundary of El sa is represented by a sequence of circular arcs. 
Each circle is characterized by its center y c l and radius r 1 ; b) For a given y, the inequalities 
(|34|) . (|35|) are satisfied within the shaded region. For the suitable i, the distance to aEl sat 
is computed as T>(ij, aEl sat ) = \\y — ay^W — ar l . 



Appendix B 

Let us discuss the computation of the derivative q~-, which enters the 
formulation of the normality rule ( 123]) . Suppose / > 0. Thus, er^ s ^ and 
the radial direction R c s := o^/\\<r^\\ is well defined. Recall that the hard- 
ening variables X^, X&, and R are to be held constant during differentiation. 
Having this in mind, we get 



D 

off II "eff 



do- ||of ff | 



ReB- (37) 



Next, taking into account that d ar ^ c ° & ( $ ^ = — (sin(arccos($))) 1 , we obtain 
lor / 



m d arccos( ^|^ ) _ 1 d(j^^ 



dcr dcr sin# dcr 

Moreover, since 9 G [0, 7r] is the angle between X d and R e g, we get 

• n -X"d — (X d : R c f[)R eS \\ 
sm = ii ii • 



(38) 
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Substituting this result into (1381) and taking into account that 
X d we get for ^ 



90 



•off 



l°"effll ll X d 



fXd : -R c fr)-R c fr| 



(39) 



Further, differentiating f l2~Tj) we obtain 
df{tr,X k ,X d ,R) 12 



do- 
where 

It follows from (JUD that 



(#0 + i?) 



9/(^(11^11,0), a) 



do- 



^3/2 



r D I 

effl 



K + R 



(cos(0),sin(0)). 



(40) 



(41) 



cffl 



dy 



r D I 

cffl 



K Q + R 



;-sin(0),cos(0)). 



(42) 

Next, note that the degree of distortion a is to be held constant as well, since 
it is a unique function of Xd. Thus, using the chain rule we get from (T40]) 



df(o,X k ,X d ,R) 
do 



^(K + R) 



df(y,a) 



'cffl 



cffl 



da 



'd f(y,a) dy\ dQ_ 
dy d6 I dcr 



Substituting fl37D, ([39]), and flU into this result we get for ^ 

df(o,X k ,X d ,R) f df(y,a) . x 
0^ = y — 7p (cos(0),sm(0)) J J Reff+ 



fdf(y,a) . 
— ' • si 



sin(0),-cos(0)) 



X r 



[X d : R c s)R c fi 



X r 



v X d : R cS )R cS \ 



■ (43) 



Here, the gradient of the non-dimensional overstress is computed by f[36|) . 
Finally, the normality vector tends to the radial direction as — > 0. Therefore 



df(o,X k ,X d ,R) 
do 



R cS for = 0. 
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